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Abstract 

Since nearly 10 years, it is known that inserting the permittivity of 
the Drude model into the Lifshitz formula for free energy causes a viola- 
tion of the third law of thermodynamics. In this paper we show that the 
standard Matsubara formulation for free energy contains a contribution 
that is non-perturbative in the relaxation parameter. We argue that the 
correct formula must have a perturbative expansion and conclude that the 
standard Matsubara formulation with the permittivity of the Drude model 
inserted is not correct. We trace the non-perturbative contribution in the 
complex frequency plane, where it shows up as a self-intersection or a bi- 
furcation of the integration path. 
PACS 03.70.-hk,11.10.Wx, 65.40.gd 



1 Introduction 

Since nearly 10 years, we observe a problem with thermodynamics in the theory 
of the Casimir effect. This problem appears when calculating the free energy T 
of the electromagnetic field using the Lifshitz formula in the presence of metallic 
bodies described by the Drude model. The same happens for dielectric bodies 
with nonzero dc conductivity. According to the third law of thermodynamics 
(Nernst heat theorem), one must expect the entropy 

s = -— (1) 

to vanish for T — i- 0. In fact, with the Drude permittivity0 

e° = 1 - (2) 
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^For 7 = it turns into the permittivity of the plasma model. 
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with plasma frequency Q and a relaxation parameter, 7(7), vanishing for T — )■ 0, 
the entropy takes a nonzero limit. 



T->0 



lim S = Si 



(3) 



which constitutes the violation of the third law of thermodynamics. This was 
found for plane parallel interfaces in |1] [2]. A similar violation was found for 
dielectrics with finite dc conductivity in [3] . Recently, a similar finding was made 
for a sphere in front of a plane [1]. 

It must be mentioned that in a number of papers, [SI El El E], nonzero 5*0 is 
not considered as a violation of thermodynamics or it is attributed to a perfectly 
symmetric crystal lattice with no impurities, which is not realized in nature [9]. 
However, these discussions are not the topic of the present paper. 

This violation constitutes not only a serious theoretical problem. Calculat- 
ing the free energy as described above, anomalous large thermal contributions at 
short separation appear first predicted in [lOj. In a very careful short-separation 
experiment with metal test bodies [HI [12] these have been shown to be excluded 
on a high confidence level. Similar large thermal corrections arising from account- 
ing dc conductivity of dielectrics were also experimentally excluded [IH [H] . In a 
large-separation experiment with centimeter-sized metal bodies [15] the predic- 
tions following from the use of the Drude permittivity were not supported. In 
line with this, it must be mentioned that in another large-separation experiment 
[16] the opposite result, supporting thermal corrections predicted by the Drude 
permittivity, was claimecl^. 

In the present paper we reconsider the foundation of the theory of the Casimir 
effect at finite temperature and we show in which place the standard approach 
with the Drude permittivity fails. We start from reminding the basic formulas 
of this approach. In application to two dielectric half- spaces with permittivity e 
separated by an empty gap with plane parallel surfaces the free energy can be 
expressed by the Lifshitz formula. 



where we have put /i = c = fee = 1- The prime on the summation sign means 
that the term with / = has to be multiplied by 1/2. In (jl]), the /-summation 
goes over the Matsubara frequencies. 



and the sum over i accounts for the two polarizations of the electromagnetic field. 
The reflection coefficients are 




(4) 



(5) 



^The method used in [TBJ was questioned in [TT] . 
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Ttm — 



(6) 



and the notation 



(7) 



is used, where fc|| is the momentum in the translational invariant directions in 
parallel to the interfaces. For e one has to insert an expression according to 
the model chosen. If e depends on the frequency u like ([2]) does, one needs to 
insert uj = i^i. This follows since the Matsubara representation (jlj) assumes a 
summation over the imaginary frequencies. 

It should be mentioned that the free energy (jl]) with the Drude model per- 
mittivity, eq. (|2]), inserted is real. This follows because on the imaginary 
frequency axis, the permittivity e° is real unlike as on the real frequency axis. 
Equation (jlj) is one of a number of various equivalent representations of the Lif- 
shitz formula. These may differ by a change of variables or by integration by 
parts. This formula was initially derived more than 50 years ago [IB]. At finite 
temperature, it uses the common representation of a quantum field theory at 
finite temperature after Matsubara. Looking at the history and on the many ap- 
plications, there seems to be no reason for any doubt in the validity of (jlj). There 
might be only a question about the permittivity that enters through the refiection 
coefficients. The basic idea is that (jl]) is the free energy of the electromagnetic 
field interacting with the dielectric bodies. This interaction is accounted for by 
the permittivities in the sense of macroscopic electrodynamics. Now, the permit- 
tivities can be taken from some theoretical model or from measurements, as well. 
Especially, the latter are well founded by a huge amount of data. In this way, 
inserting the Drude model permittivity should give a reasonable result. 

However, there is the much discussed question on what about the dissipation 
of energy inherent to the Drude model. This dissipation is due to physical pro- 
cesses like scattering and transfer of energy from the electromagnetic field into 
heat. In general, dissipation makes the considered model incomplete as long as 
the return of this heat to the electromagnetic field is not accounted for or a heat 
bath is used like in the Huttner-Barnett type of approach [12]. In physics, this 
return happens in a natural way since the heat is radiated back. 

In general, in the quantum field theory, it is common that the energy or the 
free energy acquires an imaginary part if the system is not closed and looses 
energy, by creation of particles for example. So, if not including the return 
of energy, one should expect the free energy to have an imaginary part. The 
Lifshitz formula with permittivity (|2]), interpreted as or derived from the 
vacuum energy, accounts only for the degrees of freedom of the electromagnetic 
field in the presence of the lossy dielectrics. By this setup, there is no return of 
energy included. As a consequence, one should expect the free energy to have an 
imaginary part. It is a merit of the formula (jl]) to have none. 

The question on how to include dissipation was discussed intensively in [20] [21] 
(see also the explanation in [22], section 7.3). The main idea is to include an 
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auxiliary field, which acts like a fluctuating source for the electromagnetic field 
and returns, in this way, the dissipated energy. With some restrictions on the 
parameters, eq. (jlj) was confirmed also for the case of dissipation [20l |21]. This 
discussion is generally considered as giving the correct result. 

In our discussion, we go a similar way by dividing the system into two. The 
first is the electromagnetic field in the presence of the dielectric described by the 
Drude model. The second is some mechanism that performs the return of the 
dissipated energy. Then, as an approximation, we consider the first part only. 
Of course, this procedure is quite unsatisfactory insofar as we do not have any 
information on how good this approximation is. Nevertheless we will be able to 
gain some information on representation (jlj). In dividing the system, we will have 
two contributions to the free energy. Now, both will have imaginary parts that 
must compensate in the sum. In this sense, it is natural to have free energy with 
an imaginary part even in equilibrium if considering the first part alone. 

In line with this setup, we make a point that the free energy must be pertur- 
bative in the relaxation parameter 7, at least in the sense that for 7 — )■ 0, the free 
energy turns into that of the plasma model. This follows from physical reasons in 
terms of the expectation that, for sufficiently small dissipation, the system should 
be described by the plasma model permittivity. We argue that in the underlying 
physics, there is no process that could create a nonvanishing contribution for 
vanishing dissipation. In the following, we will use the notion 'perturbative' in 
a generalized meaning by allowing for expansion terms like 7 In 7. The point is 
that all contributions beyond the plasma model are assumed to vanish for 7 — )■ 0. 

Below, we will see that the first part of the system is indeed analytic in 7. 
It is interesting to mention that the representation (jlj) is nonanalytic in 7. This 
can be inferred already from the known observation that for T — > the limiting 
value 5*0 of the entropy does not depend on the details of how fast 7 goes to 
zero. In terms of the representation 7(T) = 71T" + . . . for T — )■ used in |4j, the 
dependence on 71 drops out for a > 1. 

It should be mentioned that the above point might occur non-natural when 
thinking of dielectric bodies described by a fixed e. As observed in [23], for the 
TE mode, the limits T — )■ 00 and e — )■ 00 (ideal conductor limit) do not commute. 
Considering the ideal conductor as a zeroth order approximation, the dependence 
on 1/e is non-perturbative. However, this is not a good example. One can argue 
that this model is not physical since a dielectric with (e — 1) not decreasing for 
high frequencies cannot be realized in nature. Indeed, as soon as (e — 1) decreases 
like, for example, in the plasma model, this problem disappears. 

In the following section, we focus on the TE mode and describe how to formu- 
late the free energy in the Drude model in terms of physical frequencies and how 
to set up a perturbative expansion. In the third section, we discuss the relation 
to a representation in terms of imaginary frequencies. The conclusions are drawn 
in the last section. 

Throughout this paper we put H = c = k-g = I- 
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2 Free energy in terms of physical frequencies 



We start our discussion of free energy from first principles, not resorting to the 
Matsubara formalism. In general, the free energy T is given by the Gibbs sum, 

J^= -ThiTre"^^, (8) 

where /3 = 1/T is the inverse temperature, H is the Hamilton operator, and the 
trace goes over the corresponding Hilbert space. For bosonic excitations with 
one-particle energies resp. frequencies cuj, where J is a generic index numbering 
the excitations, this representation turns into 

^ = TX:in(2Bi„h(^)). (9) 

Also, it can be rewritten in the form 



J J 
= ^o + At-F, (10) 

where we separated the vacuum energy, Eq, which remains for vanishing tempera- 
ture. In (ITO!) . a regularization parameter s is used. We do not touch the question 
about the temperature dependence that comes into Eq through the permittivity. 
We are interested only in the second part, Aj-J-", which is the temperature- 
dependent part of the free energy. Of course, representation (|H]) is not new, see 
eq. (5.29) in [21] for example. This formula is very general and it holds for any 
bosonic system. In application to our system, we restrict our consideration to 
the excitations uj of the electromagnetic field. For these, we can write AtJ-' in 
a more specific form. 



+ 



(27r)2 [ 

J 

r^ln(l-e-^-3) A5(fc3) 
Jo ^ dks ^ 



(11) 



where, like in (jlj), A;|| is the momentum parallel to the planes and in passing from 
(fTOj) to (ITT]) . AtJ^ became the energy density per unit area. The sum over j goes 
over the discrete (for a given A;||) frequencies Uj corresponding to the waveguide 
and evanescent modes. These are all modes whose wave functions decrease in the 
directions z ±oo. The integral over accounts for the photonic modes, the 
momenta are related by 

k = ^kl + kl (12) 

and 
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is the scattering phase shift expressed in terms of the transmission coefficient, 



of the corresponding one-dimensional scattering problem for the TE polarization. 
Here, L is the width of the gap. The scattering coefficient for TM polarization 
differs by q in the parentheses to be multiplied by e. In this setup, ks is the 
momentum of the asymptotic states and 



9 = y-^^^ (15) 

is the momentum of the photons in the gap. These modes, having a real momen- 
tum ks, correspond to photons propagating in the whole space and we call them 
photonic modes. We denote their frequency by Uk^. 

The discrete frequencies Uj appear for imaginary momenta k^ = ik^ corre- 
sponding to the poles of the transmission coefficient. Here, we have 



and the Kj are solutions of the known equations, 

tan—, (symmetric solution) 



q 2 

2 



X qL . . , . ^ 

— = — cot — , (anti-symmetric solution) (17) 



with q given by 



e)kl - x2 



For real q the solutions are waveguide modes, and for imaginary g, the solutions 
are the evanescent waves. The latter appear in the plasma model and for TM po- 
larization only. The eqs. ( |T7|) are for TE polarization. That for TM polarization 
follow by multiplying the q in the denominators by e. 

Now, we turn to the frequencies of these modes. From the Maxwell equations, 
the frequencies are solutions of the equation 

ew2 = P, (19) 

where k, eq. (fT2!) . is the absolute value of the three-dimensional spatial momen- 
tum. It is essential to understand that, for a frequency-dependent permittivity, 
e(a;), eq. ( fT9i) takes the form e(ci;)w^ = k"^ and that just this equation determines 
the physical frequencies. 
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The frequencies appearing in eq. ( !TT|) are just solutions of this equation, 
whereby Uj are the frequencies of the evanescent and waveguide modes with 



momentum k = y/c^ — x| and is the frequency of the photonic modes with 

k given by eq. ( !T2l) . 

For the Drude model with permittivity (E]), we have to insert e — ?• e{uj) = 
1 — Q'^ /{u{u + i'y)) in eq. ( IT9l) . and this equation can be rewritten in the form of 
a third-order polynomial, 

{uj^-k'^){u + i-f)-con^ = 0. (20) 

It has three roots, Ua{'y, k) (a = 1, 2, 3). For not too large 7 < 2fi, the first root 
has Im(ci;i(7, k)) < and Re(ci;i(7, k)) > 0. For 7 = 0, it turns into the positive 
frequency of the plasma model. 



Wi(0,fc)|^=o = Vn2 + P. (21) 

The second root has Im(ci;2(7, fc)) < and Re(a;2(7, fc)) < 0, and it turns for 
7 = into the negative frequency of the plasma model. These two roots describe 
propagating modes that are damped according to dissipation. Their starting 
points are 

001,2(1, 0) = i (-Z7 ± V41]2 _ ^2 j . 

The third root ^3(7, k) is purely imaginary. It starts in Go's (7, 0) = 0, and it ends 
in 073(7,00) = — ^7. This root represents the over-damped mode in the Drude 
model. It has no analogy in the plasma model. We mention that these roots 
have analytic expressions, which are, however, not really helpful, and we refrain 
ourselves from displaying them. We mention only the relation 

3 

^a;,(7,fc) = 27, (22) 

i=l 

which holds for these roots. In application to free energy (fTTl) . we consider uid, k) 
as the frequency of the physical modes. 

Now, we return to the point of the perturbative expansion for small 7. From 
eq. fl20l) . it follows that the solutions Ui^, k) all have power series expansions in 
7 with a finite radius of convergence. The first few orders of these expansions are 

^3(7, A:) = -^7^^^ + 0(7^). (23) 

These are Taylor series expansions with a finite radius of convergence, i.e., these 
series do converge inside a circle, I7I < 70, in the complex plane of 7. Radius 
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Figure 1: Functions Uii^je^", k) as functions of 7 for a = 0, a = 0.27r, a = vr and 
k = 0.3. The solid lines are the real parts, and the dashed lines are the imaginary 
parts. 



7o can be found from the explicit formulas for Ui{'~f,k). Again, these formulas 
are quite inconvenient. It is easier and more instructive to consider the real and 
imaginary parts of a;j(7e*", k) as functions of 7 for several values of a (see Fig. [1]). 
For 7 < 1, the curves do not intersect. For 7 > 2, we observe bifurcations. So, we 
conclude that singularities in the complex 7-plane are all on a finite separation 
from the origin. Hence, the expansion of each Ui{'y, k) in powers of 7 has a finite 
radius of convergence. 

From the above formulas, it is seen that the perturbative expansion requires 
Q"^ + k'^ ^ 0. However, this is just the frequency in the plasma model fl2Tl) . which 
is known to be never zero. We mention that this holds not only for the photonic 
modes but also for the waveguide and evanescent modes too. At once, one can 
see from eqs. fl25]) that the coefficients decrease for A; —t- 00. This allows the 
conclusion that the expansion can be made under the sign of the integrals in 

We continue with mentioning that all functions of u entering the free energy, 
including the solutions of eq. (fT7|) and the phase shift (fT3|) . have power series 
expansions in u. This follows from their explicit form for eq.f ll3l) . or from the 
form of the defining equations (fT7|) . Now, if we insert into these the expansion of 
LJi{'y, k), we again obtain a power series expansion with a finite radius of conver- 
gence. Finally, we insert the expansion obtained into the free energy f|TT]) . Since 
the integrations converge for each order of the expansion, we conclude that AtJ^ 
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has an expansion in powers of 7 with a finite radius of convergence. 

The zeroth order of this expansion is the free energy of the plasma modeL 
The first order is proportional to ij and purely imaginary. The second order is 
proportional to (^7)^, and it is real. Higher orders behave accordingly. So, our 
conclusion in this section is that the free energy of the first of the considered 
systems goes to zero for the vanishing relaxation parameter 7. Since we expect 
the complete system to have the same property, we can conclude that the second 
part, which we do not calculate, has this property too. 

3 On the impossibility to turn to imaginary fre- 
quencies in the Drude model 

In this section, we discuss the question about turning of the momentum inte- 
gration over ks in flTT]) to the imaginary axis. The aim is to investigate how to 
pass from representation (fTTj) in terms of physical frequencies to the standard 
Matsubara formulation which is in terms of imaginary frequencies. We start 
from reminding this procedure in the plasma model, where it is possible without 
difficulty. So, we put 7 = for the moment and look on the structure of the 
representation (fTTll . We follow the procedure applied in [25]. The integral over 
^3 can be split into two integrals by writing the logarithm in ( 1T3|) as a difference 
of two. 



Now, we turn the integration path in the first contribution upwards, — )■ zx, 
and in the second contribution downwards, k — )■ — ix. We mention that there 
are no contributions from the large circles and that there are no singularities 
the paths might cross. The complex fcs-plane is shown in figure [2l The initial 
integration path goes along the real axis. It is shown by the dashed lines for the 
two contributions in eq. (I24p . On the positive imaginary axis, we have the poles 
of ^(/cs), on the negative imaginary axis, poles from t{—ks). Further, after the 



poles, there come the cuts starting from ±.i^ ~ Vl"^ — k'^y Now, the integrals 



from below the cuts have contributions from the poles that just cancel the con- 
tributions from the discrete modes, i.e., the first part in the square bracket in 
( ITTj) . After that, the integrals along the cuts remain. Here, the two logarithms 
in (jSlj) become equal - however, the other logarithms are now different since the 
frequency becomes imaginary, uj = m. the contribution turned upwards and 
uj = —iC, in the other, with 




(24) 





(25) 



For these logarithms, we note 




(26) 
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Ini(co') 



Re(w) 



Figure 2: Complex /cs-plane for contour rotation in the beginning of section 3. 



for (jj = ix and 

In (1 - e" 
for bj = —in. Function 

with 



■f^^\ = -i^ -in + ln 
' 2 



1% sm — 
2 



(27) 



(2^ 



(29) 



is a stair function appearing from passing the points ^ = irl {I integer), where 
the sinus changes sign. Obviously, is the Matsubara frequency Next, we 
observe that term in the parenthesis just cancels the vacuum energy £"0, 

which was separated in ffTOj) . So, we continue with the complete free energy J^. 
Finally, we integrate by parts and come to 



J 27r 2 [ 2 ' ^ A/^OT 



1 dfsjx) 
in dx 



In t(ix) 



(30) 



The first contribution in the curly brackets is the surface term from the {—iir)- 
contribution. Now, the derivative of the stair function gives the sum of delta 
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functions, which allow to carry out the x- integration and, as a result, the Mat- 
subara sum occurs, 

■^^^-/fiyif (31) 

This expression is known to be a version of the Matsubara sum representation of 
free energy in the plasma model. For example, it coincides up to notations with 
eq. (5.5) in [26]. Also, it coincides with (jl]) up to a contribution that does not 
depend on separation L and up to the sum over polarizations that we did not 
indicate in this section. 

In this way, we have seen for the plasma model how to pass from the rep- 
resentation of the free energy as a sum/integral over real frequencies, i.e, over 
the physical excitations, to the Matsubara representation. Now, we consider the 
question on whether this is possible for the Drude model too. We take represen- 
tation (fTT!) with the frequency ^1(7, k) as that of the physical excitations of the 
electromagnetic field, which need to be accounted for in the free energy. 

So, we start from considering solution ^1(7, k) defined in the preceding section. 
When turning k^ — )■ zx, we have to turn the complete spatial momentum 

k= ^kl + kl ^ii = i^kl - x2 (32) 

uj depends on. Therefore, we first investigate ui{'-),k) for k = ^e*" under defor- 
mation a = ... I . 

This task is a bit cumbersome, however, it does not pose any principal diffi- 
culties. For the understanding of deformation, it is necessary to consider all three 
solutions, LL!a{'y,k) (a = 1,2,3), together. These are shown (in arbitrary units) 
in figures [3] and H] as parametric plots in the complex w-plane as functions of ^ 
for several values of a. The relaxation parameter is taken to be 7 = 0.1, and the 
plasma frequency is put = 1. First, we consider a < 0.885|, which is shown 
in figure [31 Solutions ui^2 start from the points 

ui,2h, 0) = ^ (-27 ± ^402 + 72) , (33) 

whereas starts from u^l'-f, 0) = 0. These points stay fixed under the deforma- 
tion of a. For ^ — )■ 00, solutions 001^2 behave like the corresponding frequencies 
in the plasma model, i.e., 1^1^2(7, k — > 00) = ±k + . . .. The third solution takes 
imaginary values for all ^ while a = 0. It goes on the imaginary axis from to 
—27. 

Now we increase a. The curves deform and are denoted by numbers 2 to 5 
for increasing a. It is seen that the first two solutions move towards the corre- 
sponding parts of the imaginary axis, while the third solution deforms to the right 
acquiring a positive real part. At some critical value, a* ~ 0.8855 in the given 
case, the solutions Ui and touch. Beyond, for a > a*, they detach, but in a 
changed order. That means that the correspondence between starting points (for 
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Im((j) 




Re((j) 



Figure 3: Parametric plots of functions ^^(.^e"*") in the complex cu-plane as a 
function of ^ for deforming the path. Curves 1 to 5 correspond to a/(7r/2) = 
0, 0.2, 0.7, 0.85, 0.885, respectively. The curves starting from u = 0.99 — 0.05« 
represent the first solutions, Ui, that joining u = with u = —i represent the 
third solution, L03, and that in the third quadrant represent the second solution, 
U2- The value of the relaxation parameter is 7 = 1. 

^ = 0) and end points (for ^ — t- 00) changed. Hence, a turning of the integration 
path beyond a* is impossible. So, we have to conclude that within the Drude 
model it is not allowed to turn from the physical frequencies to the imaginary 
ones. Equivalently, it is not possible to come to the Matsubara representation. 

It is interesting to see what could be done about this problem. Indeed, it is 
possible to turn the contour beyond a*. For this, one needs to include from the 
beginning the contribution from the over-damped mode, u^i^^jk), in eq. ( fTT]) . 
In that case, for a > a*, we have a part of the integration path composed of 
the former ui and C03 going from 073(7,0) to ^1(7, 00 e*"). These are the curves 
in the first quadrant in figure |H For a = f the final curve of deformation goes 
just along the positive imaginary axis. This is obviously the curve used in the 
standard Matsubara representation for the Drude model (jlj). In addition, we 
have another contribution from the path joining the starting point of Ui with 
the end point of u^. These are the curves in the fourth quadrant. These are 
complex also for a = f and, if accounted for, deliver a complex contribution to 
the free energy. In the standard Matsubara representation, this contribution is 
not present. In this way, the Matsubara representation accounts for a part of the 
frequencies originating from Ui{'-f,k) and from U3{'-f,k) and disregards another 
part of these. 

In addition, we would hke to mention that solution oj3{j,k) (eddy currents) 
being inserted into free energy (fTT]) causes a problem with the convergence of 
the integrals for large momenta. First, we mention that the convergence in (fTTjl 
results from the decrease in factor ln(l — e~^'^) for u — j- 00. Now, frequency 
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Im{a») 




Re((j) 



Figure 4: Parametric plots of functions Ua{C,e~^'^) in the complex w-plane as a 
function of C, for deforming the path. Curves 1 to 4 correspond to a/(7r/2) = 

0. 886, 0.9, 0.95, 0.995, respectively. The value of the relaxation parameter is 
7 = 1. 

6^3 (7, A;) is purely imaginary. Hence, the logarithm does not decrease. One is left 
with the hope that it oscillates sufficiently fast to ensure convergence. However, 
this is also not the case since 1^3(7, A;) is bounded, < —iu3{'y,k) < 1, for real 
k. Finally, one could argue that the convergence in ( ITTl) could come from the 
decrease in the phase shift. However, the integration over k\\ remains divergent. 

In [27], the modes corresponding to solution 073(7,0;) (these are the 'eddy' 
currents), were included in addition to o;i(7,o;) in the sum over modes in ffTOj) or 
(fTTj) . However, this alone is also not sufficient to get a real free energy, resp. the 
standard Matsubara representation, what is the aim of [27j. Even for a = tt/2, 

1. e., at the final deformation of the path, from 0^3(7,0;) a part of the path is in 
the complex plane (curve 4 in Fig. 3). In addition, in ^27j, also the modes corre- 
sponding to the solutions 072(7,0;) are included. In that case, all the imaginary 
contributions compensate each other and one comes to a real free energy (already 
on the level of eq. fllOl) or ffTT]) ). However, we cannot include mode 072(7, w) be- 
cause in that case, when setting 7 = 0, we would get the plasma model twice. 
In addition, there is another problem with the modes corresponding to solutions 
072(7, w). Their negative real part would cause the sum over the eigenvalues J 
in eq. ffTOj) to diverge. Of course, this divergence does not enter the force (it 
does not depend on separation L, as can be seen in eq. Q); however, it would 
produce an unusual temperature-dependent divergence. 
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In the remaining part of this section, we will bring an argument showing that 
the Drude model inserted into the standard Matsubara formula has a contribution 
that is non perturbative in 7. For that, we consider a contour rotation like above, 
but in the opposite direction. We start from the standard Matsubara representa- 
tion (jlj) of the free energy in the Drude model. It can be represented by eq. lpTl) . 
which is the free energy for the plasma model in Matsubara representation, where 
one has to insert ([2]) at imaginary frequency, 

(34) 



The corresponding frequency condition, which gives the connection with the 
imaginary momentum perpendicular to the plane, reads 

In opposite to the third-order equation f l20|) for u, this is a second-order equation 
for X, which can be solved explicitly. 




x(0 = We^ + Y^ + ^S- (36) 



Now, we go the way back from representation fl3T|) . or from eq. (jlj), to a repre- 
sentation in terms of real frequencies. This can be done by going from eq. fl3Tj) 
backwards or by applying the Abel-Plana formula to the sum over / in ( !3T]) . The 
latter procedure was used, for example, in eqs. (11) and (13). In any case, we 
have to turn the integration path for ^ towards its imaginary axis, ^ — )■ ±iu. Let 
us follow what happens to x(^), eq. (l36j) . We show in figure j5] the deformation of 
function x(ci;e~*") in the complex x-plane for a = . . . |. Again, we take values 
7 = 0.1, k\\ = 1 and f2 = 1. For a = 0, the initial curve is along the real axis 
from 1 to 00. For increasing a, the curve goes down into the fourth quadrant. 
These are curves 1 and 2. Further increasing a, at some critical value of a, the 
curve develops a cusp and, further, a self- intersection, see curve 3. Increasing a 
further (curve 4) until the final value of 7r/2 (curve 5) the curve starts with its 
loop and, in the following course, it comes quite close (for small 7) to the real 
and imaginary axes. It is clear that, in this way, one comes to an integration over 
real frequencies u, however, one does not come to the frequencies of the Drude 
model which are the first root of eq. ( I20j) . Instead, one comes to a representation 
with integration over real u but complex = ix{—iu). This reflects the fact 
that we do not consider the complete free energy. 

In fact, this representation is not really interesting except for the appearance 
of the self-intersection and of the loop since the loop is a contribution which is 
non-perturbative in 7. We demonstrate this in figure |6l which displays the final 
path (for a = vr/2) for decreasing values of the relaxation parameter 7. The path 
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Figure 5: Parametric plot of function x(,^e *°'), eq. f l36p . in the complex Air- 
plane as a function of ^ for deforming the path. Curves 1 to 5 correspond to 
a/(7r/2) = 0.2, 0.4, 0.85, 0.95, 1, respectively. 




Re(/f) 



Figure 6: Parametric plot of function >c[—i^), eq. f l36|) . in the complex x-plane 
as a function of C, for decreasing values of 7. Curve 1 corresponds to 7 = 0.1, the 
curves 2,3,4,5 to 7 = 0.06, 0.03, 0.01, 0.001, respectively. 
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comes closer and closer to the real and imaginary axes, however, the loop remains 
and it stabilizes at a finite size. The contribution from the parts of the path that 
are along the real and imaginary axes deliver, for 7 = 0, just the free energy of 
the plasma model. Thereby, the contributions from the real axis dropped out. 
The contribution from the loop does not disappear for any finite 7, no matter 
how small, and it constitutes just the non-perturbative contribution. 

This non-perturbative contribution can also be calculated analytically using 
the methods of In fact, in [1], limit T — )■ was investigated, however, it 
turned out that limit 7 — )■ (with T fixed) can be treated by the same method. 
We delegate this calculation to the Appendix. The result is 

A^jrDrude. ^ 0) - ArJ^P^"''"" = ^^/°(0), (37) 

where function /°(0) given by eq. (55) in |1] and by eq. ( H7|) . Equation (!37|) is, 
of course, just the contribution violating the third law of thermodynamics. 



4 Conclusions 

In the foregoing sections, we considered free energy using the permittivity of the 
Drude model. We restricted our consideration to that part of the free energy, 
which results from the excitations of the electromagnetic field in the presence of 
the dielectric bodies and, in this way, we ignored the return of the dissipated 
energy. In this approach, the free energy, or more exactly, a part of it, is given 
by eq. f lTTj) and it has an imaginary part. We have described the physical modes 
of the electromagnetic field and identified that which enter ffTTl) . Next, we have 
shown that this free energy has a perturbative expansion for small relaxation 
parameter 7. Combining with the discussion in the Introduction, where we argued 
that the complete free energy must have such an expansion, we can conclude that 
the second part of the free energy, which we did not consider in this paper, also 
must have a perturbative expansion in 7. This is the main conclusion from section 
2. 

In section 3, we discussed the relation between representation f llip of the free 
energy in terms of physical frequencies and the standard representation (Lifshitz 
formula) @, which is in terms of imaginary frequencies. We started from remind- 
ing how this transition goes in the plasma model, where it is possible without 
any difficulties. Then we tried to do the same in the Drude model. Here, we 
observed a principal difficulty from the behavior of the integration path under 
rotation in the complex plane. At some angle of rotation, a kind of bifurcation 
appears and the path looses its uniqueness. From here, we conclude that it is 
impossible to pass to the Matsubara representation, at least by contour rotation. 
Further, we considered the question what happens if ignoring this problem. It 
turns out that in doing so one comes to the standard Matsubara representation 
if integrating along a part of an integration path belonging to the over-damped 



16 



mode in the Drude model and, if, at once, dropping a contribution from some 
part of the initial integration path belonging to the physical frequencies. 

We continued in section 3 with a consideration of the above procedure the 
other way round. We started from the standard Matsubara representation (jl]) 
and tried to rotate the integration path back to physical frequencies. Needless 
to say that this works fine in the plasma model. However, in the Drude model, 
under rotation, the integration path develops a self-intersection, see figure [51 As 
a result, the integration path contains a loop in the complex plane, which gives 
a finite contribution no matter how small the relaxation parameter 7 is made. 
From here, we have two conclusions. First, and not unexpected from the above 
discussion, we conclude that one cannot pass from the Matsubara representation 
to a representation in terms of physical frequencies. Second, we conclude that the 
Matsubara representation contains a non-perturbative contribution. Also, this is 
not unexpected after the discussion at the end of the Introduction. 

We mention that the procedure, described in the above paragraph, results in 
the know representation of the Lifshitz formula in terms of real frequencies (see, 
e.g., eq. (12.37) in [24J, or eq.f HOj) in the appendix) for the Drude model too. 
This is, however, different from the representation in terms of physical frequencies 
discussed in section 2. 

Finally, we consider the question whether in the perturbative approach there 
is a way to go to the Matsubara representation. In section 2, we have seen that 
the free energy of the first part of our system has a perturbative expansion for 
small 7. Thereby, we have seen that factor + fc^ appears in the denominator. 
This is seen explicitly in eq. (123|) . From this, the natural conclusion follows that 
perturbative expansion is not possible in a representation where this parameter 
is small. Fortunately, in a representation in terms of physical frequencies, this 
parameter never becomes small. Next, we mention that the plasma model is the 
zeroth order approximation for the Drude model. In the plasma model, we have 
seen in section 3 that, in passing from physical to imaginary frequencies, i.e., 
when turning the integration path, this path goes finally along the imaginary 
axis starting from ^ = 0. However, in this point, the perturbative expansion 
is not possible. From here, we draw the conclusion that in the perturbative 
approach to the Drude model, it is impossible to pass to the standard Matsubara 
formulation. This backs the above conclusion that (j4]) contains a non-perturbative 
contribution. 

It is also possible to localize the non-perturbative contribution to by ob- 
serving that one can pass from (jl]) to ffTTj) if dropping the contribution resulting 
from integration over the mentioned loop (see Fig. |6]). This contribution is com- 
plex and it has a finite, nonzero limit for 7 — >■ 0. At the end of section 3 we have 
shown the corresponding analytic expression. This is just the contribution from 
the loop in the integration path, which causes the trouble with thermodynamics. 

From the above, our conclusion is that for the Drude model one does not have 
a Matsubara representation. The commonly used one (jl]) cannot be correct since 
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it contains a non-perturbative contribution. With eq. (fTTj) . one has a part of the 
correct representation. However, one needs to add the other, still unknown, part 
that is responsible for the return of the dissipated energy to the electromagnetic 
field. 

We would like to mention two papers, [2H] and |29], where it was found that 
the complete system is perturbative in 7 in the sense we discussed above. In these 
papers, dissipation is included by coupling of the polarization fields (oscillators 
constituting the medium) to a heat bath like in the Huttner-Barnett model. The 
equations of motion for the quantized heat bath field are solved and inserted into 
the Maxwell equations. As a result, dissipation and noise polarization appear. 
Then, the energy is divided into two parts, Wi and W2, in the notations of [29] . 
corresponding to the electromagnetic field with dissipation (Wi) and the noise 
part {W2) which provides the return of energy. That part of the free energy, 
which we discuss in this paper, just corresponds to Wi. Now, in that papers, 
it is shown that in the final answer one can tend 7 to zero and arrives just on 
the original system, without dissipation. In this sense, the models considered in 
these papers are perturbative in 7. Regrettably, these papers are of restricted 
applicability. In |28] a (l+l)-dimensional model is considered, and in |29] the 
model is considered for a homogeneous dielectric medium filling the whole space. 

In the present paper, a number of questions are left unanswered. The most 
important one is about the derivation of the Lifshitz formula (jl]) in case of dissi- 
pation. In light of the forgoing discussions, we have to expect an inconsistency at 
some place in the standard derivation. It is clear this is a quite strong statement 
that needs further investigation. 

Another opened question is about a dielectric with dc conductivity, which is 
known to have similar problems with the third law of thermodynamics. Here, 
we expect a similar problem with the Matsubara representation. The point is 
that for dc conductivity, representation ffTTj) is perturbative in conductivity cr 
as can be checked easily. On the other hand, in the Matsubara representation, 
the dependence on a is clearly non perturbative. So, the prediction is that the 
transformation from physical to imaginary frequencies has a loop hole somewhere 
like that we found in this paper for the Drude model. 

Finally, we return to the discussion in the Introduction of considering the first 
part of the free energy as an approximation to the complete one. We have seen 
that this allowed us to get qualitative conclusions about the standard formulation. 
However, it is not possible to draw quantitative conclusions on how good the 
approximation is. This means we can calculate the free energy ffTTj) numerically. 
We can do that also perturbatively. The zeroth order is the plasma model, 
the first order in 7 is purely imaginary and it describes dissipation. It will be 
compensated if accounting for the second part of the system describing the return 
of energy. The second order in 7 is real, and it gives a part of the first correction 
beyond the plasma model to the complete free energy. However, the second part 
is still unknown. Moreover, there are no a'priori reasons to expect its contribution 
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to be small as compared to the first one. Hence, the only conclusion we have so 
far, is that the correction is second order in 7. Since in the experimental situations 
7 is small (the smallest of all dimensional parameters), this information might be 
helpful. 
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Appendix 

In this appendix we calculate the limit 7 — > of the free energy of the Drude 
model and prove eq. ( l37l) . We start from eq. (jlj) representing the free energy 
jrDrude ^^^q standard Matsubara representation with permittivity (E]), of the 
Drude model inserted. Technically, we follow the calculations in eqs. (4)-(16). 
Applying the Abel-Plana formula we split the free energy, 

into the vacuum energy (this is the Matsubara sum substituted by an integral) 
and the temperature-dependent part. 



00 



Function ip{C,) is given by 

/■oo 

<^(^) = / dkiikii ln(l-rV'^^) 
^0 



(40) 



for the real argument and ip{±ix) is its analytic continuation. 

As compared to [1] we changed the notations slightly. So in ^ is now k\\ 
and q in [4] is now rj. The reflection coefficient is given by 



^^Lzl (41) 



1+7/5 

with 



|2 

\- f]'^ + r] 



e + kl (42) 
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The important step for the following is to divide the integration region in the 
analytic continuation ip{ix) into two parts. The first is k\\ G [0,x] and the second 
is k\\ G [x, oo). 

In the first part, we observe that r, eq. (HTl) . can be expanded in powers of 
7. This expansion is in fact in powers (^) . However, since k\\ < x, there is no 
singularity in the integrations at least up to the first order {n=l). From the zeroth 



order (n = 0), we get just the plasma model contribution (note r] = i^x"^ — /cj^ 

here). This can be seen using a procedure similar to that in the beginning of 
section 3 and the difference is only in the variables used. 
In the second contribution. 



00 



^T-Fl^"'"'^ = — dx^^^ I {^^{ix) - ^^{-ix)) . (43) 

we made in the part of the function ^{ix), fj40|) . resulting from k\\ G [x,cxd), a 
change of variable for 77 = y^fc^"--~r^ (which is real) and denote the result by 
ip2{ix). It can be written in the form 

/■oo 

ip2{ix) = dr]r] \n{l- r^e'^'^^) , (44) 
Jo 

where now 

— ^ 7]'^ — 7] 



^+Jf^- • • . (45) 
TTW^ + + V 

The integrals in fH3|) with (jUj) inserted do converge. If we put 7 = in (jS]), we 
get formally zero since v?2(±^a;)|^^(, is real. In this way, /S.t^2^^'^'^ 1 eq. f H3|l . does 
not appear in the plasma model at all. 

The picture is different if we consider (H3ll for non-zero 7 and perform limit 
7 — )■ 0. This limit cannot be performed under the sign of the integrals directly. 
The expansion of the intergrand would start with '-y/ix and the x-integration 
would linearly diverge for a; — ?■ 0. However, following the equation that goes in |1] 
without number between eqs. (51) and (52), one can substitute the x-integration 
for a integration by a; = 7^. One comes to the representation 

-Drude ^ I 7^ ! 



=^2]^ ^^^^^ ' ^^'^'^^^ " ^^(-^70) . (46) 

From (jH]) and fHSj) . it is seen that ip2{ilC) is in fact independent on 7. Now the 
limit 7 — )■ can be performed in under the sign of the integral, and we get 

-Drude/ , ^ f d( 



A^J----(7^0) = ^^(^2(nC)-^2(-nC)), 

^ r(0), (47) 



167r2L2 

which justifies fl37|l . We mention that the rhs of eq. fj471) is, up to notations, just 
eq. (55) in [4j. 
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